


if trend ==1 
    MainDataRaw = MainData;
    MainData    = zeros(T,vars);
    for i = 1:vars
        Y = MainDataRaw(:,i);
        X = [ones(T,1), (1:T)'];
        coef = (X'*X)\(X'*Y);
        MainData(:,i) = Y - X*coef;
    end
end


obs = T - lags;    %eg, if numobs=100, then obs = 96.

% Output per worker
Y = MainData(lags+1:end, 1);            %eg, if numobs=100, then Y is row #5 through row #100.
%
dataX = zeros(obs,lags,vars);           % eg, if numobs=100 and lags=4, then each "page" of dataX is 96 x 4 
for i = 1:vars
    xvar = MainData(:,i);
    for j = 1:lags
        newx = xvar(lags+1-j:end-j);    % eg, if numobs=100 and lags=4, then first lag is row #4 through row #99.
        dataX(:,j,i) = newx;      
    end
end
%
%
% FIRST STAGE - OUTPUT PER WORKER RESIDUALS
%
% Bivariate- 
X = [dataX(:,:,1), dataX(:,:,2), ones(obs,1)];  % lags of p,u
pcoef2 = (X'*X) \ (X'*Y);
peps2 = Y - X*pcoef2;
if uselagflows ==0
    pcoef2 = repmat(pcoef2, 1,vars-2);
    peps2  = repmat(peps2, 1, vars-2);
else
    pcoef2= zeros(3*lags + 1, vars-2);
    peps2 = zeros(obs,vars-2);
    for i = 1:vars-2 
        % Loop over first stage models.  All include lags of p, u. 
        % Each loop selects one other control: f, ee, q, h, eu, rs.
        X = [dataX(:,:,1), dataX(:,:,2), dataX(:,:,i+2), ones(obs,1)];
        pcoef2(:,i) = (X'*X) \ (X'*Y);
        peps2(:,i) = Y - X*pcoef2(:,i);
    end
end
%
% Trivariate - 
X = [dataX(:,:,1), dataX(:,:,2), dataX(:,:,3), ones(obs,1)];    % lags of p,u,f
pcoef3 = (X'*X) \ (X'*Y);
peps3 = Y - X*pcoef3;
if uselagflows ==0
    pcoef3 = repmat(pcoef3, 1,vars-3);
    peps3  = repmat(peps3, 1, vars-3);
else
    pcoef3= zeros(4*lags + 1, vars-3);
    peps3 = zeros(obs,vars-3);
    for i = 1:vars-3
        % Loop over first stage models.  All include lags of p, u, and f. 
        % Each loop selects one other control: ee, q, h, eu, rs.
        X = [dataX(:,:,1), dataX(:,:,2),  dataX(:,:,3), dataX(:,:,i+3), ones(obs,1)];
        pcoef3(:,i) = (X'*X) \ (X'*Y);
        peps3(:,i) = Y - X*pcoef3(:,i);
    end
end
%
%
% 
% SECOND STAGE 
obs2 = obs - lags +1;       % eg, if numobs=100, then obs2 = 100 - 4 - 4 +1 = 93.
%
% Build new data matrix. There are 18 variables:
%   7 (u,f,ee,q,eu,h,rs)    [Drop p at this point]
% + 6 (bivariate residuals, f,ee,q,h,eu,rs)  
% + 5 (trivariate residuals,  ee,q,h,eu,rs)
%
xx    = MainData(lags:end-1, 2:end);    % eg, if numobs=100 and lags=4, then xx is rows #4 through #99.
peps  = [peps2, peps3];                 % eg, if numobs=100 and lags=4, then peps is based on rows #5 through #100.
newxx = [xx, peps];
varsx = vars - 1;
varseps2 = vars - 2;
varseps3 = vars - 3;
varstot = varsx + varseps2 + varseps3;
%
dataX = zeros(obs2,lags,varstot);          
for i = 1:varstot
    xvar = newxx(:,i);                  
    for j = 1:lags
        newx = xvar(lags+1-j:end+1-j);  % eg, if numobs=100 and lags=4, then first lag of u runs from rows #7 to #99.
                                        %     but the contemporaneous APL residual runs from rows #8 to #100.
        dataX(:,j,i) = newx;      
    end
end
dataxs = dataX(:,:,1 : varsx);  % vars are (u,f,ee,q,h,eu,rs).
dataeps= dataX(:,:,varsx+1 : varstot);
%
%
% Bivariate Regressions
newY   = MainData(2*lags:end, 2:end);   % newY is contemporaneous values of time series, (u,f,ee,q,eu,h,rs)
ucoef2 = zeros(3*lags+1, varseps2);
ycoef2 = zeros(3*lags+1, varseps2);
% NOTE - Each regression includes unemployment lags as well 
%        as lags of one other variable in dataXs (eg, f) and, 
%        correspondingly, contemporaneous and lagged values
%        of the productivity residual derived using *THAT 
%        OTHER* variable (eg, residual of p based on 
%        lags of p, u, and f). 
%
for i = 1:varseps2
    % Unemployment regression: 
    Y = newY(:,1);
    X = [dataeps(:,:,i), dataxs(:,:,1), dataxs(:,:,i+1), ones(obs2,1)];
    ucoef2(:,i) = (X'*X)\(X'*Y);
    %
    % Regression for outcome i+1:
    Y = newY(:,i+1);
    ycoef2(:,i) = (X'*X)\(X'*Y);
end
%
%
%
% Trivariate Regressions:
ucoef3 = zeros(4*lags+1, varseps3);
fcoef3 = zeros(4*lags+1, varseps3);
ycoef3 = zeros(4*lags+1, varseps3);
for i = 1:varseps3
    % Unemployment regression: 
    Y = newY(:,1);
    X = [dataeps(:,:,varseps2+i), dataxs(:,:,1), dataxs(:,:,2), dataxs(:,:,i+2), ones(obs2,1)];
    ucoef3(:,i) = (X'*X)\(X'*Y);
    %
    % Job finding regression: 
    Y = newY(:,2);
    fcoef3(:,i) = (X'*X)\(X'*Y);
    %
    % Regression for outcome i+2:
    Y = newY(:,i+2);
    ycoef3(:,i) = (X'*X)\(X'*Y);
end

    

        







